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Abstract 

The influence of molecular vibration on the Seebeck coefficient is 
studied within a simple model. Results of a scattering theory ap- 
proach are compared to those of a full self-consistent nonequilibrium 
Green's function scheme. We show, for a reasonable choice of parame- 
ters, that inelastic effects have non-negligible influence on the resulting 
Seebeck coefficient for the junction. We note that the scattering the- 
ory approach may fail both quantitatively and qualitatively. Results 
of calculation with reasonable parameters are in good agreement with 
recent measurements [P. Reddy et al.. Science 315, 1568 (2007)] 

1 Introduction 

Development of experimental techniques for constructing and exploring molec- 
ular devices has inspired extensive theoretical study of charge transport in 
molecules, with potential application in molecular electronics. [1-4] One im- 
portant issue related to stability of potential molecular devices involves heat- 
ing and heat transport in molecular junctions. This topic has attracted 
attention both experimentally [5-10] and theoretically. [11-20,22] Another 
closely related issue involves thermoelectric properties of such devices. While 
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thermo-electricity in bulk is well studied, corresponding measurements in 
molecular junctions were reported only recently. [23, 24] 

Electron-vibration interactions in the junction (leading to inelastic ef- 
fects in charge transport [21]) may cause junction heating and affect its heat 
transport properties. [22] Theoretical considerations of thermoelectric prop- 
erties so far either completely disregarded such effects (in treatments based 
on the Landauer theory [24-26]) or included it within a scattering theory 
framework. [27, 28] The latter treats the effect of vibrations as an inelastic 
electron scattering process, and changes in the non-equilibrium distributions 
of electrons and vibrations are not described in a self-consistent manner. 
It has been shown that such changes may have qualitative effects on the 
transport. [29, 30] Note also that scattering theory approaches may lead to 
erroneous predictions due to neglecting effects of contacts' Fermi seas on the 
junction electronic structure. [31] 

This paper is motivated by recent measurements of the Seebeck coeffi- 
cient in molecular junctions. [24] Considerations based on Landauer formula 
were employed there to interpret experimental data. Our goals here are: 1. 
to show the importance of vibrations for the Seebeck coefficient and 2. to 
include vibrations in a fully self-consistent way within the nonequilibrium 
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Green's function approach for calculating the Seebeck coefficient. The struc- 
ture of the paper is the following. Section [2] introduces the model, discusses 
methods used in calculations, and presents a simple (approximate) analytical 
derivation to illustrate change in the Seebeck coefficient expression (as com- 
pared to Landauer-based expression) when vibrations are taken into account. 
Section H] presents numerical results obtained in a fully self-consistent way. 
Section [5] concludes. 

2 Model 

We consider a simple resonant-level model with the electronic level |0 > cou- 
pled to two electrodes left (L) and right (R) (each a free electron reservoir at 
its own equilibrium). The electron on the resonant level (electronic energy 
Eq) is linearly coupled to a single vibrational mode (referred below as primary 
phonon) with frequency c^o- The latter is coupled to a phonon bath repre- 
sented as a set of independent harmonic oscillators (secondary phonons) . The 
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system Hamiltonian is (here and below we use ^ = 1 and e = 1) 
H = EqC^c + Skclck + (Vkcic + h.c.) 

ke{L,R} ke{L,R} 

+ uoa^a + J2 ^pb^bp + M^QaC^c + ^ U^QaQp (1) 

/3 f3 

where (c) are creation (destruction) operators for electrons on the bridge 
level, (cfc) are corresponding operators for electronic states in the con- 
tacts, a"!" (a) are creation (destruction) operators for the primary phonon, 
and fejj (bp) are the corresponding operators for phonon states in the thermal 
(phonon) bath. Qa and are phonon displacement operators 

Qa = a + a^ Qf3 = b(3 + bl (2) 

The energy parameters Ma and Uf^ correspond to the vibronic and the vi- 
brational coupling respectively. Eq.([T]) is often used as a generic model for 
describing effects of vibrational motion on electronic conduction in molecular 
junctions. [21] 

After a small polaron (canonical or Lang-Firsov) transformation [32,33] 
the Hamiltonian ([1]) takes the form (for details see Ref. [29]) 

H = eoC^c+ ekclck+ (^Vfcc|,cX„ + h.c.) 

ke{L,R} ke{L,R} 

+ uod^a + ^ cJ^b^b/s + ^ UpQaQp (3) 
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where 



M2 

So = So -A (4) 

LUq 



Xa = exp 



iXaPa] Xa=— (5) 



A is the electron level shift due to coupling to the primary phonon and Xa is 
primary phonon shift generator. Pa — —i{d — a^) is the phonon momentum 
operator. 

The mathematical quantity of interest is the single electron Green func- 
tion (GF) on the Keldysh contour 

G{ti, T2) = -i<n c(ti) c\t2) >h (6) 
Following Ref. [29], we approximate it by 



= G,{n,T2)}C{n,T2) (7) 



where Gc{ti,T2) is the pure electronic GF while IC{ti,T2) corresponds to 
the Franck-Condon factor. In Ref. [29] we have developed a self-consistent 
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scheme for evaluating this functions, leading to the coupled set of equations 

/C(n, T2) = exp {a^ [thDp^pSn, T2)- < Pi >] } (8) 
Dp^pSry) = Df^,S^y) (9) 
+ jdri jdT2 Dflp^{T, Ti) np„p^(ri, T2) Dp^p^{t2, t') 

G,{t,t') = G^^\t,t') (10) 
+ Yl /^^i / (r, n) ^,,K{n, r2) G,{r2, r') 

where Dp^p^{t,t') = —i < TcPair) Pair') > is the phonon GF, Dp}p^ and 
are the phonon and electron Green functions when the two sub-systems 
are uncoupled (M^ = 0). The functions IIp^p^ and T,c^k in Eqs. ([9]) and (fTOi) 
are given by 

Up^P^in,T2) = J2\Up\'Dp,P,irur2)-tXl Yl l^^l' 

X [figk{T2, Ti) Gc(ri, Ta) /C(ri, rg) + (n ^ T2)] (11) 
) = Y \^k\'9k{ri, r2) < T,Xa{r2)Xi{n) > (12) 

These functions play here a role similar to self-energies in standard many- 
particle theory. Here K = L,R and Qk is the free electron Green function for 
state k in the contacts. For details of derivation see Ref . [29] . A self-consistent 
solution scheme implies solving Eqs. fl51)- ffT^ iteratively until convergence. 



As a convergence parameter we used population of the level no =< CqCq >. 
When no for subsequent steps of the iterative cycle differed by less than a 
predefined tolerance (taken in the calculations below to be 10"^), convergence 
was assumed to be achieved. 

Once the electron GF ([7]) is obtained, its lesser and greater projections 
are used to get steady-state current through the junction [34, 35] 

Ik = (Ik) = J^ [S< (i?) G>{E) - T>{E) G<{E)] (13) 

at interface K = L,R. Here 

E<(^) = ifK{E)TK{E) (14) 
S>(E) = -i[l-fKiE)]TKiE) (15) 

with fxiE) = [exp(/5(£' — ^k)) + 1]^^ the Fermi distribution in the contact 
K = L,R and 

rK{E) = 2nY,\ykME-ek) (16) 

keK 

The Seebeck coefficient is defined by 

SH) ^ ^ (17) 

where V{I) is the voltage bias that yields current / at AT = 0, while AT(/) 
is the temperature difference between contacts that yields the same current 
at V = 0. The linear regime corresponds to the / — > limit of Eg. (1171) . 
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Below we present calculations of the Seebeck coefficient using different 
levels of approximations. In particular we compare results of a simple scat- 
tering theory-like approach and a full self-consistent calculation based on the 
procedure described above. The simple approach is essentially a first step of 
the full self-consistent iterative solution with the additional assumption of no 
coupling to the thermal bath for the molecular vibration (f/^ in Eq.([n])). 

3 Transport coefficients 

Before presenting results of numerical calculations, we point out how trans- 
port coefficients are introduced. In the Landauer regime of transport (electron- 
phonon interaction disregarded, both carriers scatter ballistically) , the elec- 
tric and thermal fluxes, I and J, are given by [22,25,36] 





(18) 




where Ep is a common Fermi energy in the absence of bias. 



%i,o{E) = Tt[Tl{E) G\E) Vr{E) G"(E)] 



(19) 



rp,,o(^) = TA^dE) D'{E) QniE) D\E)] 
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are the electron and phonon transmission coefficients in tlie absence of electron- 
plionon coupling, and where 

nf(u;) = 27rJ2\Uf3\'^^{i^-^f^) K^L,R (20) 

is the broadening of the molecular vibration due to coupling to its thermal 
environment. In the linear response regime the currents are linear in the 
applied driving forces - the bias V and the temperature difference AT 

I^G-V + L-AT 

(21) 

J^R-V + F-AT 
Here G and F are the electrical and thermal conductions, respectively, and 
L is known as the thermoelectric coefficient. The coefficients are given by 
[25, 36] 

G = -— y dE [-(3f\E)]%i^o{E) (22) 

I I p+oo rp Tp 

L^-^j ^dE [-f3f{E)]%,,o{E)—^ (23) 
lei /■+°° 

^=-^y dE[-(3f{E)]%i^o{E){E-EF)^L-T (24) 



-oo 

f JP Z7'_^2 



1 r+oo 

F= — / dE[-(3f'{E)]%i,o{E) 



{E - Ef) 



T 

— oo 

foo , ,2 



1 fco 2 

+ 2^1 ^'^M^'(^)]V,o(a;)y (25) 
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where f'{E) is derivative of Fermi-Dirac distribution, N'{uj) is derivative of 
Bose-Einstein distribution, T is temperature [jS = 1/T), and Ep is the Fermi 
energy in the leads. Note the existence of Onsager relation, L-T = R, between 
the cross coefficients. Note also that the coefficient F in fl25p contains two 
contributions, one corresponding to energy transfer by electrons, the other 
- by phonons. For discussion of the additive from of F and the relative 
importance of these contributions see Ref. [22]. The Seebeck coefficient is 
given in terms of these transport coefficients by 

S = ^ (26) 

Below we focus on these two coefficients - G and L - only. Making the 

approximation [—f3f'{E)] ^ 5{E — Ep) in fl22|) . and utilizing the Sommerfield 

expansion [37] in ( !23l) . Eq.( !26l) leads to 

^ _ n^k^T d\n%i,o{E) 
^ ~ dE 

which is Eq.(4) of Ref. [25]. 

4 Calculation of the Seebeck coefficient 

As discussed in Section [21 the simplest calculation that takes into account the 
electron- vibration interaction term (the Ma term in Eq.([T])) corresponds to 
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inelastic scattering of the transmitted electron from the phonon at the given 
initial temperature. Within the self-consistent scheme presented in Section [2] 
this result is obtained after the first iteration step, with the influence of the 
(free) vibration on the electronic GF is taken into account but not vice versa. 
For this reason the coupling f/g to the thermal bath can be disregarded. For 
the model this calculation yields [39] 
I I +00 

^ = ^ y. 4(2AViVo(iVo + l))/„^(2AVA^o(iVo + 1)) 

n,m=— 00 

/+00 
dEZi^oiE) (28) 
-00 

X {hiE + nuo)[l - fR{E - muo)] - fn^E + muo)[l - UiE - nuo)]} 

where A'^o = [e^'^" — and /„ is the modifled Bessel function of order n. [38] 
For Ma = Eg. (128!) reduces back to (fT8|) . Linearization in the bias potential 
V = Vl — Vr and in the temperature difference AT = Tl — Tr leads to the 
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phonon-renormalized transport coefficients 

2 +°° 

G = -^ V 4(2AVA^o(iVo + l))/™(2AVA^o(iVo + 1)) 
nil ^ — ^ 

n,m=— oo 

/ + 00 
dE%i^o{E) (29) 
-oo 

X {[-/3/'(^ - m^o)]/(^ + nujo) + [-f3f'{E + muo)][l - f{E - nuo)]} 
L = -^ 4(2AV^o(iVo + l))/™(2AV^o(iVo + 1)) 

n,m=— oo 

^ g/3(n+H-o/2-2A2(2JVo+l) / dE%i^o{E) 



X |[-/3r(E-m^o)]/(i^ + ^^o) ^ ""^^ (30) 
+ [-mE + mu,m - fiE - nc.o)] ^ + "^;°~^" 

Using Eqs. (I29l) and (l30l) . the Seebeck coefficient is calculated from Eg. (126 



To this end we ffist calculate the currents I{V, AT = 0) and I{V = 0, AT) 
as functions of V and AT. The inverted functions V{I, AT = 0) and 
AT {I, V = 0) are then used in fl2Bl) to yield S{I) (expressed below as 
S(y) with V = V{I,AT = 0)). In the calculations presented below we 
have used symmetric bias and temperature differences across the junction: 
fiL,R = Ef± V/2, Tl,r = T± at/2. 

Figure [1^ shows the Seebeck coefficient 5* as a function of the bias po- 
tential, calculated at T = 300 K using the energetic parameters Ep = 0, 
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Figure 1: (Color online) Seebeck coefficient vs. applied bias. Shown are 

results of full self-consistent calculation (solid line, red), scattering theory 

approach based on model (128!) (dashed line, blue), and elastic scattering case 

(dotted line, black) calculated with (a) 'standard' set of parameters (see text 
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for parameters), (b) higher Eq — Ep gap, and (c) weaker electron-phonon 
coupling Ma- 



eo = 0.2 eV, = = 0.005 eV, ujq = 0.05 eV, M„ = 0.1 eV, and 
Qf = QP^ = 0.005 eV. The last is a wide-band approximation for molecu- 
lar vibration broadening (|20|) due to coupling to thermal baths (for detailed 
discussion see Ref. [40]). Figures [TJd and c show similar results for a higher 
£0 — Ep gap and a weaker electron-phonon coupling Ma, respectively. To 
show the effect of electron-phonon coupling on the Seebeck coefficient we 
compare the results obtained from the full self-consistent calculation and the 
scattering theory approximation to the elastic Ma = limit. 

In the following figures we consider the inelastic effects only within the 
scattering theory approximation (which requires a far smaller numerical ef- 
fort). The dependence of S on the energy gap Eq — Ep is shown in Fig. [21 and 
its variation as function of ujq is displayed in Fig. [31 In these figures V is kept 
at the value 0.05 V and all unvaried parameters are the same as in Fig. [H 
Finally, in Fig. [H we compare self-consistent results (some of them already 
shown in Fig. [ID obtained for different choices of electron-phonon coupling 
Ma and vibrational broadening VL^^. 

The following observations can be made regarding these results: 

1. In contrast to inelastic tunneling features usually seen in the second 
derivative (Pl/dV"^ of the current- voltage characteristic near \eV\ = 
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Figure 2: (Color online) Seebeck coefficient vs. energy position of the molec- 
ular level for the model of Eq. (!28|) . Shown are results with (dashed line, blue) 
and without (dotted line, black) electron-phonon interaction. Calculation is 
done at V = 0.05 V. Other parameters are the same as in Fig. [1] 
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Figure 3: (Color online) Seebeck coefficient vs. vibrational mode frequency 
for the model of Eq. (!28|) . Shown are result with (dashed line, blue) and 
without (dotted line, black) electron-phonon interaction. Calculation is done 
at = 0.05 V. Other parameters are the same as in Fig. [1] 
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Figure 4: (Color online) Seebeck coefficient vs. applied bias for model de- 
scribed in Section [21 Solid line (red) and dash-dotted line (magenta) are 
identical to solid lines (red) in Figs. [1^ and c respectively. Also shown result 
for Q^j'^ = = 0.002 eV (dash-double dotted line, black). Elastic case 
(dotted line, black) is shown for comparison. 
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hujQ no such threshold behavior is seen in Fig. [TJ This lack of threshold 
behavior in the inelastic contribution to the Seebeck coefficient results 
the fact that thermoelectric conduction is assosiated with the tails of 
the lead Fermi-Dirac distributions, and these tails wash away any ther- 
shold structure. 

2. Inelastic contributions can have a substantial effect on the Seebeck co- 
efficient and its voltage dependence, however the assessment of these 
contributions is sensitive to the approximation used and cannot gener- 
ally be based on the scattering theory level of calculation. Indeed, as 
seen in Figure [T] S{V) behavior may change qualitatively upon going 
from the scattering to the self-consistent calculation. 

3. Focusing on the self-consistent results. Figs. [T] and H] show that the 
inelastic effect on the Seebeck coefficient can be positive or negative, 
depending on the other energetic parameters in the system. A change 
of sign as a function of Eq is seen also in the scattering theory results of 
Fig. [21 The existence of a similar crossover in the self-consistent results 
can be inferred by comparing Figures [T^ and b. 

4. A smaller value of the electro n-phonon coupling Ma naturally leads 
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to a smaller difference between self-consistent and scattering theory 
calculations (compare Figs. [1^ and c). Still the coupling chosen in 
Fig. [TJ: is strong enough to give appreciable difference between the 
inelastic and elastic results. 

5. In Fig. [2] the Seebeck coefficient is seen to go through a maximum as 
a function of the gap Eq — Ep, with inelastic contribution affecting the 
position and height of the observed peak. The behavior seen in Fig. [2] 
can be rationalized by noting that for, say, Ti > Tji electrons with 
energies E > Ep contribute most to left-to-right current, while those 
with E < Ep dominate right-to- left current. This gives no thermo- 
electric current when Eq = Ep, hence as Eo ^ Ep one needs a higher 
temperature difference in order to compensate for the same bias. As a 
result, the Seebeck coefficient goes down at — Ep — ^ 0. On the other 
hand when Eq — Ep ^ F, ksT the two contributions cancel each other, 
and hence Seebeck coefficient drops down once more. 

6. The dependence of 5* on uq (Fig. [3] demonstrates it within a scattering 
theory level calculation) is in line with the expectation that S should 
attain its classical limit as the vibration becomes more rigid. 
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7. We have found (not shown) that effects on S of varying the electronic 
(r) or vibrational {Q) widths show a similar trend as varying the gap 
Eq — Ep or the vibrational frequency tuo, respectively. 

5 Conclusion 

We studied the influence of molecular vibration (inelastic effects) on the See- 
beck coefficient for molecular junction transport using a simple model of one 
molecular level (representing the participating molecular state) coupled to 
two contacts and to one molecular vibration. Two approaches to the model 
were considered: a simplified scattering model represented by Eg. (128!) and the 
full self-consistent approach described in Section [21 The simplified approach 
ignores the mutual infiuence of electronic and vibrational subsystems. Note, 
the structure of the expression for current, Eq. fl28|) is just a difference between 
two scattering fiuxes (left-to-right minus right-to- left). Results of simplified 
model calculation are compared to the full self-consistent approach where 
both mutual electron-vibration infiuence and vibration coupling to thermal 
baths are taken into account (for detailed description of the approach see 
Ref. [29]). We show that inelastic effects have non-negligible influence on 
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the resulting Seebeck coefficient for the junction, for a reasonable choice of 
parameters. Electron-vibration interaction can either increase or decrease 
Seebeck coefficient depending on the physical situation. We study the de- 
pendence of this influence on different parameters of the model (applied bias, 
gap between molecular level and Fermi energy, strengths of coupling between 
molecule and contacts, between tunneling electron and molecular vibration, 
between molecular vibration coupling and thermal baths, and frequency of 
the vibration). Comparing results of the two approaches, we show that scat- 
tering theory based approach may fail both quantitatively and qualitatively. 
Experimental data presented in Ref. [24] does not give conclusive evidence 
for the relative importance of inelastic processes. More extensive measure- 
ments showing dependence of the Seebeck coefficient on junction parameters 
are needed in order to make definite conclusion. In particular, isotopic ef- 
fects should influence vibration-related part of the Seebeck coefficient (see 
Figs. [3]and|4]). Change in electron-phonon coupling would also reveal the 
inelastic part of the Seebeck coefficient (seeHj). Finally, results of our model 
calculations done with a reasonable set of parameters provide the Seebeck 
coefficient to be of the order of ~ 10~^ V/K. Results reported in Ref. [24] 
yield S ~ 10~^ V/K. Since molecules used in the experiment [24] are charac- 
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terized by relatively big gap Eq — Ep, our results are in good agreement with 
the measured data. Indeed, for eq — Ef ~ 1 eV calculated Seebeck coefficient 
(see Fig. [2]) becomes of the experimentally observed order of magnitude. 
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